Molecular Dynamics

The following tutorial runs two molecular dynamics (MD) simulations with BigDFT and analyses its outputs. We investigated the dynamics of a water environment with NVE ensamble at DFT level and the motion of a 38 atoms Lennard-Jones cluster with fixed NVT.

from BigDFT import Logfiles as lf
from BigDFT import Calculators as calc
from BigDFT import Inputfiles as I
import matplotlib.pyplot as plt

A water environment

To study and reproduce bulk water, we simulate the motion of four water molecules placed in a periodic cubic cell. The cubic unit cell side has been chosen to adjust the effective density of water to its ambient density.

n = 4 # number of water molecules
m = 18.01528 # [g/mol] molar mass
NA = 6.02214076e23 # [1/mol] Avogadro number
rho = 997e-27 # [g/AA^3] water density at room temperature

volume = ( n * m / NA ) / rho
box_size = volume ** (1.0/3.0)
print('Box size: {} [AA]'.format(box_size))

Box size: 4.932703172202558 [AA]

from io import StringIO
input_file = StringIO("""HETATM    1  H   HOH     0       2.924   3.219   1.550  1.00  0.00           H
HETATM    2  H   HOH     0       4.098   2.445   0.948  1.00  0.00           H
HETATM    3  O   HOH     0       3.850   2.984   1.714  1.00  0.00           O
HETATM    4  H   HOH     1       1.432   4.761   0.988  1.00  0.00           H
HETATM    5  H   HOH     1       0.066   4.168   0.637  1.00  0.00           H
HETATM    6  O   HOH     1       0.884   3.972   1.119  1.00  0.00           O
HETATM    7  H   HOH     2       2.280   0.928   1.171  1.00  0.00           H
HETATM    8  H   HOH     2       0.907   1.567   1.386  1.00  0.00           H
HETATM    9  O   HOH     2       1.335   0.779   1.019  1.00  0.00           O
HETATM   10  H   HOH     3       4.694   4.067   3.659  1.00  0.00           H
HETATM   11  H   HOH     3       3.377   4.728   4.071  1.00  0.00           H
HETATM   12  O   HOH     3       3.945   4.569   3.303  1.00  0.00           O  """)

from BigDFT.Atoms import AU_to_A
from BigDFT.IO import read_pdb
from BigDFT.UnitCells import UnitCell

sys = read_pdb(input_file)
sys.cell = UnitCell([box_size, box_size, box_size], units="angstroem")

from BigDFT.Visualization import InlineVisualizer
viz = InlineVisualizer(400,300)
viz.display_system(sys)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Geometry optimization

To avoid large forces (and associated large spatial displacements) at the first MD step, we relaxed the initial geometry of the atomistic system.

Soft norm-conserving pseudo-potentials including non-linear core corrections along with the Perdew-Burke-Ernzerhof functional were used to describe the core electrons and exchange-correlation both for the geometry optimization and the MD trajectory. The wavelet basis functions were distributed on an adaptive uniform mesh with a resolution of \(h_{\text{grid}} := h_x = h_y = h_z = 0.40\) Bohr. The \(h_{\text{grid}}\) of a single KS SCF calculation has to be chosen to guarantee a well converged energy as done for an usual DFT single point calculation.

water_inp = I.Inputfile()
water_inp.set_hgrid(0.4)
water_inp.set_rmult([8.0, 10.0])
water_inp.set_xc("PBE")
water_inp.set_symmetry(False)
water_inp.optimize_geometry(method="SBFGS", nsteps=500, frac_fluct=3.0, forcemax = 1e-4, betax=1.0)
# water_inp['geopt'] = {'method': 'SBFGS', 'ncount_cluster_x': 500, 'frac_fluct': 3.0, 'forcemax': 1.E-4, 'betax': 1.0E0}
study = calc.SystemCalculator(skip=True)

Initialize a Calculator with OMP_NUM_THREADS=2 and command mpirun -np 2 /Users/wddawson/Documents/CEA/binaries/bds/install/bin/bigdft

water_geopt = study.run(posinp=sys.get_posinp(), name="water-geopt",
                                       input=water_inp) #Run the code with the name scheme geopt

Creating the yaml input file "./water-geopt.yaml"
Executing command:  mpirun -np 2 /Users/wddawson/Documents/CEA/binaries/bds/install/bin/bigdft -n water-geopt -s Yes
Found 71 different runs

water_geopt.geopt_plot()

We can extract the steps of the geometry optimization process, and then visualize.

from copy import deepcopy
optsys = []

for inst in water_geopt:
    optsys.append(deepcopy(sys))
    optsys[-1].update_positions_from_dict(inst.log["Atomic structure"])

from BigDFT.Visualization import InlineVisualizer
viz = InlineVisualizer(400,300)
viz.display_system(*optsys)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol